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ABSTRACT 

The binary star system y Cephei is unusual in that it harbours a stable giant planet around 
the larger star at a distance only about a tenth of that of the stellar separation. Numerical 
simulations are carried out into the stability of test particles in the system. This provides 
possible locations for additional planets and asteroids. To this end, the region interior to the 
planet is investigated in detail and found to permit structured belts of particles. The region 
between the planet and the secondary star however shows almost no stability. The existence 
of an Edgeworth-Kuiper belt analogue is found to be a possibility beyond 65 au from the 
barycentre of the system, although it shows almost no structural features. Finally, the region 
around the secondary star is studied for the first time. Here, a zone of stability is seen out to 
1.5 au for a range of inclinations. In addition, a ten Jupiter-mass planet is shown to remain 
stable about this smaller star, with the habitability and observational properties of such an 
object being discussed. 

Key words: stars: individual: y Cephei (HR 8974) - planetary systems - binaries: general - 
methods: A^-body simulations 



1 INTRODUCTION 

7 Cephei is one of the closest separation binary systems that con- 
tains a planet. The primary and secondary stars are separated by 
18.5 au. The planet (designated y Cep Ab) has a minimum mass 
of 1.7 Mi and follows an eccentric orbit about the larger star 
alone. The observat ional parameters of y Cephei, as determined by 
lHatzes et al] 120031) . are summarised in TableQ 

Although the data are reasonably reliable, there is some uncer- 
tainty regarding the masses of the components. The inclination to 
the line of sight of a spectroscopic binary is usually indeterminate. 
This means that the mass of the primary is generally reckoned from 
spectral type. Then, the mass function (e.g., Smart 1977) permits 
a minimum mass to be assigned to the secondary, albeit crudely 
(Griffin, Carquillat & Ginestet 2002). 

For 7 Cephei, the mass of the prima ry is » 1.57M Q from 
photospheric modelling (Fuhrmann 2004). The most recent de- 
termi nation from stellar evolution models is x \.lM a lAffer et al. 
l2005t) . The mass of the secondary is given bv lDvorak et alJ 12003) 
as ss 0.4M Q . However, assuming the original value of the primary's 
mass of 1.57 M Q , a minimum mass can be derived as « 0.34 M Q 
from the velocity amplitude fitted bv lHatzes et alJ 1200 3l) . 

There have been a number of studies of the dynamics of the y 
Ceph ei system to d ate. Foremost is the numerical investigation of 
lDvoraketalU 2003). which does use an earlier, and slightly differ- 
ent, set of orbital parameters (see Table 0. They used Burlisch- 
Stoer and Fast Liapunov Indicator methods to show that the y 
Cephei system is stable over Myr time-scales. They carried out test 
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particle integrations as a guide to the possible existence of further 
planets. The results show a stable inner region between 0.5 and 0.8 
au and then an a dditional stable zone at low inclination around 1 
au. lDvoraketal1l2003l) noted that this coincides with the 3: 1 mean 
motion resonance. This resonance is stable in the y Cephei system, 
but is unstable in the Solar system. [Dvorak et"aH 120031) conclude 
that Earth-mass planets (up to 90 M E ) can exist in this region which 
is fortunately j ust on the edge of t he habitable zone 1 . An extension 
to this work by Dvorak 'et alJl2004l) showed that the planet's eccen- 
tricity is less important than that of the two stars for the stability of 
a second planet. HashiahiDour] 120051) also carried out numerical 
studies of the system's stability using a Burlisch-Stoer integrator, 
for a range of possible configurations of the planet's and binary's 
semimajor axes, eccentricities and inclinations, confirming and ex- 
tending the results of lDvorak et alji2003l) . 

Here, we use both numerical simulations and analytic calcula- 
tions to study zones of stability as possible locations of additional 
planetary companions to either star and asteroid or Edgeworth- 
Kuiper belt analogues. Edgeworth-Kuiper belts are of particular in- 
terest in binary systems, as they may be being observed (indirectly) 
in exosystems such as r Ceti and rj Corvi (Greaves et al. 2004; 
Wyatt et al. 2005). The results fall into three categories: possible 
planetary companions and asteroid belts centred on the primary star 
(§2), planetary companions around the secondary star (§3) and fi- 
nally possible Edgeworth-Kuiper belt about both stars (§4). 



The habitable zone has boundaries which depend on assumptions re- 
garding stellar luminosity and effective temperature, as well as the climate 
model adopted for the hypothetical habitable planet. 
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Name 


y Cep A 


y Cep B 


y Cep Ab 


Class 


K1IV sub-giant star 


M dwarf star 


Planet 


Mass 


1.59 ±0.12 M 


0.4 M e 


1.7 + 0.4 Mj 


Period (days) 




20750.6579 ± 1568.6 


905.574 ± 3.08 


Semimajor axis (au) 




18.5 +1.1 


2.13 ±0.05 


Eccentricity 




0.361 ± 0.023 


0.12 ±0.05 


Longitude of periastron (° ) 




158.76 ± 1.2 


49.6 ± 25.6 


Time of periastron passage (JD) 




2448429.03 + 27.0 


2453121.925 ±66.9 



Table 1. Best fit orbital parameters for the y Cephei system from Hatzes et al. r2003). Mass of star B is from lDvorak et alj 120031) . 



Name 


y Cep A 


yCepB 


y Cep Ab 


Class 


KIIV sub-giant star 


M dwarf star 


Planet 


Mass 


1.6M 


0.4M o 


1.76Mj 


Period (days) 




25567.5 


902.2 


Semimajor axis (an) 




21.36 


2.15 


Eccentricity 




0.44 


0.209 



Table 2. Orbital parameters for the y Cephei system used bv lDvorak et all 
120031) . 

2 PLANETS AND ASTEROID BELTS AROUND THE 
PRIMARY 

2.1 Algorithm 

For all the simulations, we use the parameters of the y Cephei sys- 
tem given in Table Q Here and elsewhere in the paper, the equa- 
tions of motion are integrated using a conservative Burlisch-Stoer 
metho d provided in the MERCURY software package (Chambers 
1999). Although not as fast as sympletic methods, this was chosen 
because of its ability to provide close encounter data and handle 
highly eccentric objects. The two stars and planet are simulated as 
point masses for gravitational interactions. Any additional objects 
are taken as massless test particles to decrease integration times. 
Test particles are removed from the simulations when they collide 
with the primary star, or pass an ejection distance of the order of 
several hundred au. A collision with the primary means that the 
test particle has a position that lies within the body of the primary, 
as judged by its stellar radius of 0.02 au (Hatzes et al. 2003). Close 
encounters are allowed to occur, and are defined as taking place 
whenever a test particle enters within one Hill radius Rh of the sec- 
ondary star or the planet, defined as 

/ m \ 1/3 

where m is the mass of the secondary or planet and M is the mass 
of the primary. This works out as » 8.1 au for the secondary and 
« 0.15 au for the planet. 

To maintain accuracy, the variation in the system's total energy 
and angular momentum is monitored throughout each simulation. 
Using this to constrain the initial timest ep to 1 day and the tolerance 
in the Burlisch-Stoer algorithm JPress et aljl 999) to 10~ 12 leads to 
an overall fractional change in the system's energy AE/E of about 
10~ 8 over a 100 Myr period. This can therefore be considered the 
maximum time the system can be accurately followed. All the sim- 
ulations presented here are typically 1 Myr in time-scale, for which 
AE/E ss 10"" or better. 

It is straightforward to show that the y Cephei system has long 
term stability. A 100 Myr integration shows no major variation 
in the orbits. Regular short period variations do occur, for exam- 
ple, a slight oscillation of the planet's semimajor axis over time- 
scales equal to both its orbital period and the binary's orbital pe- 




• Test particles Plane of the test particles 

Figure 1. Two possible configurations for test particles. Top: the planet 
remains in the plane of the binary and the test particles are inclined relative 
to this. In other words, the inclination of the planet (^b is the same as the 
inclination of the secondary ;b and both are zero. Bottom: the planet and 
test particles have the same inclination relative to the plane of the binary. In 
other words, iAb is equal to the inclination of the test particles i tp . 

riod. An additional secular variation is seen over a period of about 
5500 years, evident in the eccentricity and longitude of the planet 
only. This secular period is in good agreement with the results of 
quadrupole theory (see eq. (36) of Ford et al. (2000)). 

2.2 Test Particles Interior to the Orbit of the Planet 

Figure shows two possible configurations of test particles in the 
system considered here. In the first, the test particles are inclined 
to the common plane of the binary and planet. In the second, both 
the test particles and planet are coplanar, yet inclined relative to 
the plane of the binary. The binary is always assumed to be viewed 
edge-on (that is, it lies in a plane perpendicular to the plane of the 
sky). If the planet is at an inclination i Ab relative to this, its mass 
must be increased accordingly by dividing by sin i where i = 90° - 

; Ab- 

With the planet and binary coplanar, we begin by investi- 
gating the stability of test particles in the region interior to the 
known planet. A grid of particles with semimajor axis from 0.5 
to 1.85 au and inclination from 0° to 50° with resolution 0.05 au 
and 5° respectively is integrated for 1 Myr. Thirty-six particles are 
started at each grid point with varying initial longitudes of peri- 
centre (a) = 0°, 60°, . . . 300°) and longitudes of ascending node 
(CI = 0°, 60°, . . . 300°). The orbits are initially circular. The sta- 
bility is then determined by computing the mean survival time (t s ) 
in Myrs at each grid point averaged over the 36 test particles. The 
results are shown in Figure [2] and can be compared to figure 2 of 
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Figure 2. Stability map for test particles interior to the orbit of the planet 
in the coplanar case i^b = 0. Here, stability is indicated by the test particle's 
survival time averaged over longitudes and normalised to 1 Myr. Darker 
colours indicate more stable regions, whilst lighter colours show less sta- 
ble regions. Over-plotted in purple are the nominal locations of the mean 
motion resonances with the planet Ab, up to the 5:4 case. 



iDvorak et alj 120031) . Note th at our stability index is slightly dif- 
ferent to the criteria used bv IDvorak et all l2003h . who removed 
test particles after they become orbit crossing. This may miss the 
occasional test particle that is stable, for example, if it lies in a 
Trojan-like orbit. We only remove test particles if they collide with 
the central star or are ejected from the system. 

The map shows test particles with semimajor axes less than 
« 1.4 au are stable. However, there is a strip of instability between 
roughly 0.8 and 1.0 au, creating an island of stability at low inclina- 
tions between 1.0 and 1.3 au. Some of the structure in the map can 
be clearly related to the positions of the mean motion resonances 
(MMRs) with the planet (indicated in Fig.0- The 4:1, 3:1, 5:2 and 
2: 1 resonances seem to mark transitions from stability to instabil- 
ity. For example, the 5:2 MMR divides the island at a 1.15 au. The 
lack of effect of some of the higher order MMRs, such as the 5:1 
case, is probably due to their comparative weakness and narrow- 
ness. The lack of stability beyond as 1.4 au is readily explained. 
The gravitational reach of the planet as a multiple of the Hill radius 
(Jones, Underwood & Sleep 2005) places the limit on stability at 
1.31 au, a good match with the results here. Note that in calculating 
this limit the maximum eccentricity of the planet obtained during 
the simulation has been used. Many of the test particles in the high 
inclination region of Figure [2] show ev idence of Ko zai cycles, as 
illustrated by the orbits in Figure^ The lKozaH ^62) instability is 
well-known from studies of high inclination comets and asteroids 
in the Solar system. It sets in at inclinations greater than a critical 
value of ; cl it = asinV0~4 =s 39.23°. During Kozai cycles, the ec- 
centricity and inclination vary so as to maintain approximate con- 





Figure 3. The variation of eccentricity and inclination with time for a test 
particle undergoing Kozai cycles. The starting elements are a = 0.5 au and 
i = 50°. 



I ll-: [a = 0°) 



3:1 MMR (o = 90°) 




Figure 4. The variation of eccentricity with time for test particles at the 3:1 
MMR for starting a) = 0° and 90°. 



stancy of the integral of motion / K = ^1 - e 2 p cos «, p , where e lp is 
the test particle's eccentricity. As the semimajor axis increases, the 
amplitude of eccentricity and inclination librations becomes larger, 
thus accounting for the increased instability evident in this region 
of Figure|2| 

As already seen, the MMRs with the planet are important 
in shaping the regions of stability. However, one major feature 
unexplained by this is the instability strip between roughly 0.8 
and 1.0 au. Although at zero inclination the edges are marked 
by the 4:1 and 3:1 MMRs, the instability strip shows a pro- 
nounced evolution with inclination which is suggestive of a sec- 
ular resonar^ejnsteadTh^ classical Laplace-Lagrange linear the- 
ory (Murray & Dermott 2000), although derived for low eccentric- 
ity and inclination regimes around a dominant central mass, can be 
applied to give a first approximation of the locations of these reso- 
nances. A secular resonance for a test particle occurs when its pre- 
cession rate has exactly the same magnitude as an eigenfrequency 
of the system. The eigenfrequencies are easily calculable for the 
three body system made up of the two stars and planet and are the 
eigenvalues of the 2x2 matrices A and B respectively, which have 
components 

. 1 ™-k .,(1). . 

Ait = aab\Aa), 

" 'AM + trij 3/2 



Bjj 



-Hi — aab?L (a), 

J 4M + mi 3/2V ' 



-A, 



(2) 



b\-; 2 (a) 
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Here, n, is the mean motion of object j (1 represents the planet 
and 2 represents the secondary), m ; and aj are the mass and semi- 
major axis of object j, M is the mass of the primary (the cen- 
tral object), a = a l /a 2 and b 3 ) 2 (a) and by 2 (a) are Laplace coef- 
ficients. Using the values for the masses and semimajor axes given 
in Table along with n\ = 145.2° yr 1 , n 2 = 6.337° yr 1 and 
a = 0. 1 15 1 gives the Laplace coefficients as b^) 2 (a) = 0.3542° yr 1 

a nd b?Ua) = 0.05089° yr' employing the MATHEMATICA routines 
of Mu rray & Dermottl (2000). Calculating the matrices and solving 
for the eigenfrequencies gives 

gi = 0.04338° yr 1 , 

g 2 = 0.00005211° yr 1 , (3) 
/ = -0.04343° yr 1 , 

where g\ and gi are the eigenvalues of A and / is the degenerate 
eigenvalue of B. Note that the g t and / eigenfrequencies have about 
the same magnitude, whilst g 2 is almost zero due to the large mass 
ratio between the planet and secondary star. The precession rate of 
the test particle is given by 

1 lm\ _ m m 2 _ ni \ 

A,p = "I [~M aiSlb Vi ( - ai) + j^^/jfe)] (4) 

where n is the particle's mean motion. For the region interior to the 
planet atj = atj = a/a,, where a is the test particles semimajor axis. 
Plotting A tp as a function of a from 0.5 to 1.85 au shows a resonant 
location where the gi and / eigenfrequencies intersect the curve 
at « 0.8 au. This supports the idea that the location of the inner 
edge of the instability strip on the map coincides with a secular 
resonance. 

Comparing our results with those o flDvoraketaTi r2003). it is 
easy to see that the broad trends are similar, despite slightly differ- 
ing orbital eleme nts. The main stable regions are slightly closer to 
the star in lDvorak et alJ 120031) . This may be due to a higher ec- 
centricity of both the planet and the secondary, which means that 
they approach the ce ntral star more close ly, reducing separations 
with the test particles. iDvorak et alJ <2003l) find that the 3: 1 MMR 
is stable, in contrast to asteroids in the Solar System in the same 
resonance with Jupiter. Here, we find that the resonance is unsta- 
ble, with a particle following a fairly steady evolution until switch- 
ing to a mode where its eccentricity is rapidly driven to unity on a 
time-scale of 10 kyrs, as shown in Figure|4] 

The case where the planet is also inclined (i Ab t 0) has not 
been previously studied. This is a more likely configuration for a 
system that has formed in a common protoplanetary disc. Here, 
we investigate this case by using the same grid of test particles as 
before, but with the planet sharing the same inclination as the test 
particles and with Qi* = 0. This means that, to reproduce the same 
radial velocity dataset, the mass of the planet must be increased by 
dividing by sin i, where i = 90° - i A t- 

The results are displayed in Figure|5] As compared to the ear- 
lier case of Figure|5| the unstable region has expanded, especially at 
high inclinations. This is partly caused by the change in the extent 
of the gravitational reach of the planet, as shown by the dotted curve 
in Figure|5| Although this does scale with the increasing planetary 
mass, the sharp change at 40° inclination is due to the planet be- 
coming subject to Kozai cycles. The large increase in eccentricity 
here means that the planet's periastron is much closer to the star, 
and hence its gravitational influence is larger. At 50°, the planet is 
unstable, colliding with the central star after as 0.5 Myr. 

At first sight, it may seem that the rest of the increased in- 
stability evident in Figure|5|is caused by the increased mass of the 
planet. However, experiments show that this is not so. For example, 
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Figure 5. As Figure[2] but here the planet and test particles are all inclined 
to the plane of the binary. Over-plotted in purple are the nominal locations 
of the mean motion resonances with the planet Ab, up to the 5:4 case. The 
dotted curve shows the gravitational reach of the planet, calculated using 
the method of Jones et al. (2005). 



an integration of the 30° case with the planet inclined but at min- 
imum mass (1.7Mj) shows very little difference to Figure [5] This 
suggests that the cause lies in the relative inclination of test parti- 
cle and planet to the plane of the binary. In the case of Figure [2] 
the amplitude of libration of ! tp and e tp of test particles is generally 
modest for all cases below i crit , the critical value for the Kozai in- 
stability. For test particles in Figure [5] the amplitude is no longer 
small and becomes larger with increasing inclination, rapidly driv- 
ing particles into the regime where the Kozai instability is effective. 
This is understandable as in the former case, the forces due to the 
masses in the system are always directed towards the plane of the 
binary, whereas in the latter case this is not true. As the inclina- 
tion increases, the misalignment between the forces due to the stars 
and the force due to the planet increases, and so the amplitude of 
libration increases. 



2.3 Test Particles Interior to the Orbit of the Binary 

Holman & Wiegert (1999) have already studied the stability of test 
particles in binary systems. These may orbit either a single star or 
both stars. Here, we study test particles around the primary star. 
For this case, Holman & Wiegert (1999) introduced the notion of 
a critical semimajor axis a crit , which is the largest circular orbit 
around the primary for which a ring of test particles survives for 
at least 10 4 binary periods (« 600 kyrs in the case of y Cephei). 
Using their eq (1), we find a clit = 4.0 ± 0.6 au. In other words, our 
expectation is that all test particles starting out at semimajor axes 
greater than 4.0 au will be rapidly swept out. The existence of the 
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Figure 6. The ultimate fate of test particles interior to the orbit of the binary, 
with inclinations 0° , 30° , 1 50° and 1 80° . For clarity, only test particles in 
and above the plane of the binary and planet are displayed. The following 
abbreviations are used in the key: UE = unstable, ejected from the system; 
UA = unstable, collided with star A; UE(Ab) = unstable, ejected from the 
system after a close encounter with Ab; UA(Ab) = unstable, collided with A 
after a close encounter with Ab; UE(B) = unstable, ejected from the system 
after a close encounter with B; UA(B) = unstable, collided with star A after 
a close encounter with B; S = stable; S(Ab) = stable, but suffered a close 
encounter with Ab; S(B) = stable, but suffered a close encounter with B. 



comparatively large and eccentric planet y Cephei Ab will cause 
further de-stabilization. 

To investigate this, simulations of test particles in the region 
from 0.5 to 20.0 au are carried out. They have initially zero ec- 
centricity and are set up on a grid with resolution 0.5 au in start- 
ing semimajor axis. The range of inclinations is restricted from 
0° to 30° for the prograde case, and from 150° to 180° for the 
retrograde case, both in steps of 10°. Seventy-two particles are 
started at each grid point with varying initial longitudes of peri- 
centre (a> = 0°, 30°, . . . 330°) and longitudes of ascending node 
(£2 = 0°, 60°, . . . 300°). The orbits of the test particles are followed 
for 1 Myr, and any close encounters are recorded. Although this is 
a limited range of inclinations, it is expected that those not investi- 
gated are largely unstable. This receives confirmation from explo- 
ration integrations in the case of 70° inclination, for which very few 
test particles survive throughout the entire region. 

Figures |6| and show the ultimate fates of the test particles 
within 5 au. They differ in that the planet is in the plane of the 
binary and test particles are inclined in Figure|6| whilst both planet 
and test particles are similarly inclined in Figure Four panels 
showing the results of selected simulations are displayed in each 
case, corresponding to prograde with / tp = 0° and 30°, retrograde 
with ( tp = 180° and 150°. 

In the inner regions, test particles close to the planet are swept 
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Figure 7. As Figure |6| but the test particles and the planet are similarly 
inclined. 



out on a precession time-scale («; 5.5 kyrs). Within as 3 au, pro- 
grade test particles are either ejected or collide with the central star 
after a close encounter with the planet. This is evident from the 
particles coloured orange and red in the upper panels of both Fig- 
ures |6| and The retrograde case is different. The lower panels of 
Figure|6|show swathes of stable particles coloured either brown or 
green, according to whether they reside within the Hill sphere of 
the secondary or not. Note that the secondary's periastron is at 1 1.8 
au, so particles out to 3.7 au are within its Hill sphere at some point. 
In the case i tp = 1 80° , the retrograde stability zone extends out to 
~ 7 au. As the inclination decreases (i tp — > 150°), the stability zone 
shrinks to the annuli between 0.5 to 1.0 au and 3.0 to 5.0 au. The 
lower panels of Figure Q show the case when both test particles 
and planet are retrograde. In the case i lp = 180°, the large stability 
region seen previously has almost completely disappeared. There 
are only a few remaining test particles that survive the 1 Myr in- 
tegration. The similarity of the two right-hand panels of this figure 
shows that the inner region's evolution is almost entirely controlled 
by the planet. 

The numbers of test particles surviving after 1 Myr at each 
semimajor axis are given in Table [3] Shown in Figure |H| are the 
mean survival times plotted against semimajor axis for a variety 
of inclinations. The mean survival time (t s ) can be computed by 
averaging over the results for the differing longitudes of ascending 
node and pericentre at fixed semimajor axis and inclination. The 
averaging is over 12 test particles for i tp = 0° and 180° and 72 for all 
other cases. Test particles that survive to the end of the integrations 
are included with a survival time of 1 Myr, so that the computed 
mean survival time is a lower limit in these cases. Figure [8] shows 
that there is a region of enhanced stability with (t s ) as 100 kyrs - 
for all the studied inclinations - centred at a 3.5 au, just beyond the 
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Figure 8. The mean survival time of test particles as a function of starting semimajor axis for particles interior to the orbit of the binary. In the left-hand 
panels, the starting inclination of the test particles is 0° (red), 10° (yellow), 20° (green) and 30° (blue); in the right-hand panels, it is 180° (red), 170° (yellow), 
160° (green) and 150° (blue). The top panels are inclined test particles only, the bottom panels are similarly inclined test particles and planet. The starting 
semimajor axis of the secondary star and the planet are shown as vertical dashed lines. 



Table 3. The number of test particles interior to the orbit of the binary at each starting semimajor axis and inclination that survive for 1 Myrs. There are 72 
test particles initially for all cases except i = 0° and (' = 180°, which have 12. 



Starting Semimajor Axis [in au] 

Inclination 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 12.0 12.5 13.0 13.5 14.0 
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Figure 9. The mean survival times (black line) and fate of individual parti- 
cles (coloured points) for test particles around star B. The results are shown 
for inclinations 0° and 30°. The fate of each particle is given by the key: 
UE = unstable, ejected from the system; UB = unstable, collided with star 
B; S = stable; UA = unstable, collided with star A; UE(Ab) = unstable, 
ejected from the system after a close encounter with Ab; UA(Ab) = unsta- 
ble, collided with star A after a close encounter with Ab; S(Ab)= stable, 
but suffered a close encounter with Ab. Note that, in the retrograde cases, 
almost all particles are stable for the length of the integration (1 Myr). 

gravitational reach of the planet but within the critical semimajor 
axis and just within the region affected by the secondary. Most test 
particles here (and also beyond this region) suffer a close encounter 
with the secondary before ejection or collision with the primary. 
The effect of the secondary star becomes increasingly important 
with increasing semimajor axis and (/ s ) falls to a 1 kyr. There is, 
for the prograde cases, a region of enhanced stability around 12 au 
for some inclinations, clearly visible in Figure [8] This is due to a 
few particles surviving here for the full length of the integration. 

As stated in the introduction, the mass of star B is uncertain. 
By altering this parameter and rerunning some of the simulations 
the importance of the uncertainty can be seen. For the coplanar con- 
figuration of test particles described in this section, changing the 
mass of star B by ~ 20% results in almost no difference in global 
statistical properties, such as average survival times. This would in- 
dicate that the system's dynamics are not significantly affected by 
the uncertainty in the secondaries mass. 



3 PLANETS AND ASTEROID BELTS AROUND THE 
SECONDARY 

Here, we investigate whether planets could exist around star B. 
This possibility has not been investigated before for this system. 
The critical semimajor axis for stability, as defined by Holman & 
Wiegert (1999), is 1.9 ± 1.0 au. To investigate further, test particles 
are set up from 0.5 to 2.5 au in steps of 0.05 au for inclinations 0° 
to 30° and 150° to 180° in steps of 10°. The particles are on ini- 
tially circular orbits again, and spaced in longitude of perihelion 
and ascending node by 60° as before. The results of two of the pro- 
grade cases are shown in Figure Prograde test particles are not 
stable beyond 1.5 au. In the case of 10°, 20° and 30° inclinations, 
the unstable test particles generally survive at least ten times longer 
than those in the 0° case. The unstable test particles within about 
2 au all collide with the primary, whilst those beyond this point 
have a range of fates. Holman & Wiegert's critical semimajor axis 




0.5 1.0 1.5 2.0 2.5 0.5 1.0 1.5 2.0 2.5 

nitiol semimojor axis (au) Initial semimajor axis (au) 



Figure 10. The average maximum eccentricity and maximum change in 
semimajor axis for test particles around star B. This is calculated by aver- 
aging over longitude at each semimajor axis. The inclinations shown in the 
top panels are: 0° (red), 10° (yellow), 20° (green) and 30° (blue). The in- 
clinations shown in the bottom panels are 180° (red), 170° (yellow), 160° 
(green) and 150° (blue). Also shown as dashed vertical lines are the semi- 
major axis within which all particles remain stable. For the 150° case the 
one unstable semimajor axis is shaded in grey. 




Figure 11. The radial velocity curve (left) and residuals (right) for the 
primary star when there is an additional planetary companion to star B. 
The additional planet has a mass of lOMj . The radial velocity curve shows 
the original system in black and the new two planet system in red. The 
residuals are calculated by subtracting these two curves. When there is no 
planet around star B the star has an increased mass of Mb + Meb, where 
Mrjb is the mass of the additional planet to be added to the later simulation. 

does not match up as well as in Section l2~3l but still agrees within 
the rather large uncertainty. In the retrograde cases, almost all the 
test particles are stable, with the exception of 7 particles at 2.45 au 
in the 150° case. Integrating more distant particles shows that the 
retrograde stability reaches out to about 3.5 au. 

The stability of the test particles can be further investigated 
by plotting the evolution of their eccentricity and semimajor axis, 
as shown in Fieure fTol For all the simulations, there is not much 
change in the inclination of the test particles. The variation in ec- 
centricity and semimajor axis of stable particles increases as the 
initial distance from star B increases, and is similar for all the pro- 
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grade cases. However, the variation is much larger for the (more 
stable) retrograde cases. When i tp = 150°, the variations are no 
longer smoothly increasing and show abrupt jumps, indicating that 
this case may not be long term stable. 

As test particles can remain stable around the secondary, this 
raises the question of whether a massive planet could also persist 
here. So, it is interesting to consider whether this could be de- 
tectable in the radial velocity curve of the primary star. To investi- 
gate this, the simple case of a 10 Jupiter-mass planet in a coplanar, 
circular orbit with initial semimajor axis 0.5 au and initial longitude 
0° is integrated. The orbital elements of the secondary are adjusted 
so that the centre of mass of it and its putative planet orbits the 
primary with the elements shown in Table 1 for star B alone. As 
expected from the test particle results, the planet remains stable for 
the full Myr length of the integration and shows very little variation 
in its orbit. However, the 1 Myr integrations here may overestimate 
the extent of the stability zone for planets, as instabilities can appear 
even after 100 Myr of apparent stability (Jones, Sleep & Chambers 
2001). To calculate the radial velocity curve, the system is assumed 
to be at zero inclination relative to the line of sight (i = 90°). Fig- 
ure ll ll shows the radial velocity curve of the primary, together with 
the residuals with respect to the case with no extra planet. There is 
no detectable signal with the period of the extra planet. The small 
variations shown, which have the period of the binary, amount to 
25 ms 1 over 1000 yr timespans, would be undetectable. 



4 EDGEWORTH-KUIPER BELT ANALOGUES 

The final region remaining to be studied is that exterior to the bi- 
nary. This may host particles analogous to the Edgeworth-Kuiper 
belt in our own Solar system. There is observational evidence 
for the existence of extrasolar Edgeworth-Kuiper belts from in- 
frared imaging of dusty discs around other stars (e.g., Wyatt 2003, 
Greaves et al. 2004), so the longevity of such a structure around y 
Cephei is worth investigation. 

Numerical integrations of test particles around binary systems 
suggest that they will remain stable in this region (e.g., Harrington 
1977). In addition to the criterion in Section l2~31 Holman & Wiegert 
(1999) also provide an empirical rule-of-fhumb to determine pro- 
grade test particle stability exterior to the binary. Here, the critical 
semimajor axis a crit is that beyond which almost all test particles 
survive for at least 10 4 binary periods. Using eq. (3) of their paper, 
we find that a cr j t = 64 ± 2 au. However, Holman & Wiegert (1999) 
caution that there can sometimes be islands of instability beyond 
the critical radius associated with n:i MMRs. There is an older 
criterion due to Harrington (1977), who also considers the case of 
retrograde test particles. Harrington's equation suggests that pro- 
grade test particles are stable for a > 56 au and retrograde stable 
for a > 45 au. By stability, Harrington means that the particles 
show bounded motion with no secular or large periodic variations 
in their elements over his - by nowadays standards - very small 
integration timespans. Using the same method as in Sections 12.21 
and !2.3l test particles are set up with inclinations to 30° and 150° 
to 180° in the (barycentric) region from 20 to 150 au in steps of 

5 au. The mean survival times are shown in Figure [121 The pro- 
grade test particles exhibit a sharp cut-off, with those beyond 65 
au being stable, independent of the starting inclination. This figure 
also shows that the retrograde test particles survive to the end of 
the integration for starting semimajor axes beyond 40 au for 180° 
inclination, 45 au for 170° and 160° inclinations and 60 au for 150° 
inclination. The range of initial conditions for which retrograde par- 
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Figure 12. Average survival time of test particles exterior to the binary 
as a function of barycentric initial semimajor axis. In the left hand panel 
the starting inclinations are: 0° (red), 10° (yellow), 20° (green) and 30° 
(blue). In the right hand panel the starting inclinations are: 180° (red), 170° 
(yellow), 160° (green) and 150° (blue). 





20 40 60 80 100 120 140 

Initial barycentric semimajor axis (au) 



20 40 60 80 100 120 140 

Initial barycentric semimajor axis (au) 







150° 






160° 






170° 


- 

1 . 


1 


ISO" 



31 10 








% S 






160° 

170" 


E 

■ 1 6 






iao" 


CT> 4 








| 2 












i 





20 40 60 SO 100 120 140 

Initial barycentric semimajor axis (au) 



20 40 60 80 100 120 140 

Initial barycentric semimajor axis (au) 



Figure 13. The average maximum eccentricity and maximum change in 
semimajor axis of each test particle in the region exterior to the binary, as 
a function of starting semimajor axis. The inclinations shown are as for 
Figure [T2l and the averages calculated as for Figure fTol The boundaries 
between the unstable and stable regions are shown as dashed grey lines. An 
isolated unstable region is shown as a shaded area for the 150° case, as in 
Figure fTol 



tides survive is larger than for prograde (e.g., Harrington 1977). 
For all inclinations, unstable particles from about 20 to 40 au are 
removed quickly and generally are ejected after a close encounter 
with the secondary. Unstable particles from about 40 to 65 au are 
less rapidly removed and tend to collide with star A or be ejected 
from the system, without experiencing a prior close encounter. We 
see that both Harrington's (1977) and Holman & Wiegert's (1999) 
stability criteria seem to give a reasonable description of the results 
of our simulations. 

The average maximum eccentricity and change in semimajor 
axis is shown in Figure [13*1 as a function of initial semimajor axis 
for all eight inclinations. All the test particles show a small secular 
variation in their orbital elements. 
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Figure 14. The average survival times (black line) and fate of individual 
particles (coloured points) for test particles in the region exterior to the bi- 
nary. This is for an inclination of 0° and with a resolution of 1 au in semima- 
jor axis. Only the particles out to 80 au are shown, as the simulation is stable 
beyond here. The locations of the n: 1 MMRs are shown as dashed grey lines 
and labelled above the plot. The colour coded fates of test particles are as in 
Figurelol 



The Edgeworth-Kuiper belt in our Solar system shows some 
structure due to locations of MMRs with the giant planets (e.g., 
Luu & Jewitt 2002). However, the results in Figure fT2l do not in- 
dicate any such resonant features. This is most likely due to the 
coarse semimajor axis grid employed. To investigate this, the copla- 
nar case was re-run, but now with a spacing in semimajor axis of 
1 au. The mean survival times and fates of individual particles for 
this simulation are shown in Figure fT?! The first location where all 
test particles survive is now 64 au, agreeing exactly with Holman 
& Wiegert's value of a cr ; t . There is then an unstable band from 67 
to 70 au that seems to match up to the 7: 1 MMR. The plot gives no 
evidence for any other resonant effects beyond this location. This 
is understandable since the locations of the other major resonances 
are outside the region of Hill stability (for example, the 3:2 reso- 
nance that is important in our own Solar system), leaving only those 
of very high order to affect stable test particles. The distinct differ- 
ence in fates of test particles is obvious in this plot, with the blue 
coloured points, indicating close encounters with star B, not occur- 
ring beyond about 45 au. At this point, the test particles are yellow 
or grey, indicating either ejection or collision with star A without 
undergoing any close encounter. 



5 CONCLUSIONS 

In this paper, we have carried out a suite of test particle integra- 
tions for the y Cephei system, as a guide to possible locations for 
additional planets. For test particles in the plane of the binary and 
planet, there are three zones of stability for 1 Myr timespans at 
least. These are [1] the region interior to the planet, which is stable 
within the bands 1.2 and 1.3 au, 1.05 to 1.15 au, and 0.75 to at least 
0.5 au, [2] the region around star B from 1.5 au into at least 0.5 au 



and [3] the region around both stars extending out from about 65 
au. These can be used to constrain possible locations of additional 
planetary companions. 

The region interior to the planet has a complicated structure. 
Low order resonances, such as the 4: 1, 3: 1 and 5:2, mark transitions 
between stable and unstable regimes. There is a secular resonance 
located at 0.8 au from the primary that also plays a role in the dy- 
namics of the test particles. The results for this region match up 
well with the previous work of Dvorak et al. (2003) despite the 
differences in the parameters of the system and the methods used. 
This implies that the slight improvements in the observations of the 
system have not significantly altered its dynamical characteristics. 
One difference though is that, unlike Dvorak et al. (2003), we find 
that test particles close to the 3:1 resonance are unstable, just as for 
the asteroid belt in the Solar System. In agreement with Dvorak et 
al. (2003), we find that small planets interior to y Cephei Ab are 
long-lived. Another as yet unconsidered possibility seen from the 
results here is the existence of an asteroid belt interior to the giant 
planet. However, the stability seen for test particles in this region 
disappears when y Cephei Ab is inclined in the same plane as the 
test particles. Now, the test particles are rapidly driven to inclina- 
tions which are subject to the Kozai instability. 

The region interior to the binary and exterior to the planet is 
unpromising. It seems unlikely that any prograde planet or asteroid 
can remain stable between y Cephei Ab and the binary, although 
retrograde test particles in this region are more long-lived. The re- 
gion around the secondary is stable out to 1.5 au for test particles. 
We have shown that planets up to 10 Jupiter-masses in circular or- 
bits at 0.5 au can survive for at least 1 Myr. This seems to be a 
promising place for additional planets to reside, in a similar man- 
ner to satellites about a planet around a single star. 

In the final region, that exterior to the binary, test particles are 
also stable. Although planets would be able to reside out here, the 
existence of an Edgeworth-Kuiper belt structure is also a possibil- 
ity. Once again the retrograde particles are more stable than those 
that are prograde. Whilst in every region studied this is true, such 
objects are perhaps unlikely, if the origin of the system was a com- 
mon disk. If particles are captured from elsewhere, then this may 
become a possibility. Retrograde asteroids and comets certainly ex- 
ist in the Solar system. 

Since the survival of additional planets in the system has been 
shown to be a possibility, it is interesting to consider their habitabil- 
ity. There are a number of estimates of the habitable zone around 
star A in the literature. For example, Dvorak et al. (2003) place it 
at 1.0 to 2.2 au, in which case habitable planets could exist at the 
very edge of the zone. However, Haghighipour (2005) places the 
habitable zone at 3.1 to 3.8 au, while Jones et al. (2005) place it at 
2.07 to 4.17 au. These different locations reflect differences in the 
criteria for the location of the habitable zone or the assumed stellar 
luminosity and effective temperature. The zone 2.07 to 4. 17 au is 
unstable for all except retrograde test particles. So, this would per- 
mit the possibility of a habitable planet only if retrograde. The de- 
tection of any such small planet is challenging, given the small (of 
the order of a few ms~') radial velocity signatures that they cause. 
In addition, a retrograde planet appears in radial velocity curves 
merely as a prograde one with a 180° phase difference. Although 
retrograde objects - especially planets - are thought to be unlikely, 
it has been shown that giant planets in binary systems can end up 
with retrograde orbits after close encounters within a planetary sys- 
tem (Marzari et al. 2005). 

A more interesting possibility is that of a habitable planet 
around star B, since it has been shown that planets can survive 
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here. Theories of planetary formation do not preclude this possi- 
bility (e.g., Armitage, Clarke & Palla 2003). The habitable zone 
for an M dwarf extends out to about 0.3 au (Kasting, Whitmire & 
Reynolds 1993), which is within the stable zone found here. The 
large primary star nearby might also act as a 'shield' from comets 
and asteroids similar to Jupiter for the Earth. As seen in the results 
in Sections [5] and [4] very few test particles collide with the sec- 
ondary star once the region interior to the binary has been rapidly 
cleared. Edgeworth-Kuiper belt objects that are perturbed into the 
inner regions of the system are also likely to suffer encounters with 
the primary, thus leaving the secondary and its environs largely un- 
scathed. 

The detection of such a habitable planet around star B has, 
unfortunately, been shown to be virtually impossible from the ra- 
dial velocity signature of star A alone. Should the secondary be 
resolved this would change, and any companion giant planet would 
be easily detectable. The region around star B, none the less, re- 
mains as a promising place for the existence of a habitable planet 
in this system. 
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